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I. INTRODUCTION 



In this article we present a method to compute the scattering states of holes in spherical bands 
in the strong spin-orbit coupling regime. More precisely, we calculate scattering phase shifts and 
amplitudes of holes induced by defects in a semiconductor crystal. We follow a previous work 
done on this topic by Ralph [H. I. Ralph, Philips Res. Rept. 32 160 (1977)] to account for the 
p-wave nature and the coupling of valence band states. We extend Ralph's analysis to incorporate 
finite-range potentials in the scattering problem. We find that the variable phase method provides 
a very convenient framework for our purposes and show in detail how scattering amplitudes and 
phase shifts are obtained. The Green's matrix of the Schrodinger equation, the Lippmann-Schwinger 
equation and the Born approximation are also discussed. Examples are provided to illustrate our 
calculations with Yukawa type potentials. 

PACS numbers: 72.10.-d, 72.10.Fk 
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The valence bands of the semiconductors of the column IV (Si, Ge) and of some III-V (GaAs, InP) as well as II- VI 
(CdTe, ZnSe) compounds have a similar structure. Because of the complexity of this structure, calculation of the drift 
t-H , mobility fi^ and the mobility of Hall /in in these materials, when they are of the p-type, is a very difficult problem. 
Simplifying hypothesis usually consist of neglecting the split-off band and assuming parabolic band dispersions, Bj(k), 
for the light (i = £), and heavy (i = h) hole bands. Even in these conditions, calculations remain complicated since 
one has to take into account the p-wave nature of the holes wavefunctions, as well as the various possible transitions 
between the two subbands. 

In the semi-classical model, the evaluation of the transport coefficients requires the knowledge of the non-equilibrium 
\Q ' hole occupation functions /»(k) (i — £,h), solutions of the Boltzmann equation. The transition rates between the 
states \i, k) and \ j, k'), characterizingthe various hole scattering processes (phonons, defects ■ • ■ ) are included in the 
definition of the collision integral 0,12. In the present work, we focus on the problem of elastic scattering by a point 
, defect, modeled by a spherically symetric perturbation potential. 

In a weak field regime, the Boltzmann equation can be linearized and the scattering by the defects are characterized 
by four relaxation times, , that correspond to the four possible transitions within (intra) and between (inter) the 
two bands i and j. In the Born approximation, the various T}j are calculated introducing the overlap factors Gy (k, k') 
El 0,1^ that characterize the p-w&ve nature of the hole wavefunctions. For the Yukawa potential, frequently used to 
model ionized defects, expressions of the times are well known 0,0- I n the low temperature regime, where the 
scattering by defects become dominant compared to the other scattering mechanisms, the thermal energy of the holes 
• i-h , is small compared to the perturbation induced by the defects in the crystal, and the Born approximation ceases to be 
^ ■ valid. 

To overcome this problem, the phase shift method which yields an essentially exact solution to scattering problems, 
can be used. In the case of carriers that belong to the same band (no coupling) of s-wave nature, the expression of the 
relaxation time r as function of the phase shifts 5i of each partial wave with angular momentum I, is well known jg]. 
Meyer and Bartoli have used this result to perform approximate calculations of carrier mobility limited by scattering 
on ionized impurities in p-doped Si and GaAs 9]. In their model they neglected the interband transitions and ignored 
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the p-wave nature of the holes wavefunctions. They considered a Yukawa potential and gave approximate analytic 
forms allowing direct calculation of mobilities without the need to compute the phase shifts. 

Kim and Majerfeld have calculated phase shifts to compute the scattering states of holes on a single ionized 
impurity jlOj . Their model has nonetheless two shortcomings: they assume a weak perturbation and their phase 
shifts calculations are performed using a Hamiltonian different from the well established one usually found in many 
articles reporting on the defects in p-type semiconductors 0, El E3> 0] • 

A more rigorous treatment of the scattering of holes based on the phase shift method was done by Ralph [HI- 
Calculations were done with a spherical Hamiltonian in the strong spin-orbit coupling regime. For each of the two 
regular independent solutions of the Schrodinger equation, the radial functions associated to each scattering partial 
wave were completely characterized by two amplitudes, a^, and two phase shifts, 5 ± , in the asymptotic region 
(r — > oo); and are called Ralph's parameters in the present article. In this model the contribution of each 
partial wave to the various scattering times thus required the calculation of eight parameters: four scattering 
amplitudes and four scattering phase shifts. The example treated by Ralph was only an academic one since a zero- 
range potential of the Dirac-type was used, for which it is possible to obtain simple analytic expressions of the r,j . 
Ralph's parameters had never been computed for more realistic potentials, such as a Yukawa potential, much used 
for the study of electron and hole scattering by ionized defects. 

In this article we show how to make use of the variable phase method |r| to obtain Ralph's parameters for finite 
range potentials. The variable phase method has been mostly developped and used for collisions and scattering 
problems in atomic and nuclear physics, see e.g. 11^. ilTL flSL ll9L l20t l2lL 12^. l23| . whereas there are comparatively few 
works where it is used in solid state physics, e.g. |2J,|23. To the best of our knowledge this method had never been 
applied to scattering problems with coupling in solid state physics. In its simplest version (no coupling) the variable 
phase approach consists of deriving a first-order differential equation satisfied by the phase shifts 8i (r) of the partial 
waves with angular momentum I induced by a potential truncated at distance r, and solve it imposing <5;(0) = as 
an initial condition. The method is also useful to calculate the amplitudes ai(r) of the partial waves. For the mixed 
valence-band holes the present work consists of tranforming a 2 x 2 second order differential system satisfied by the 
two radial wavefunctions into a 4x4 first order differential system whose solutions are the two amplitudes a^{r) and 
the two phase shifts 5^{r). 

The article is organised as follows: in Section II, the spherical model Hamiltonian on which we base all the subsequent 
calculations is introduced. In Section III we establish the framework of our calculations. We show how we derive the 
generalized coupled first-order differential equations satisfied by the scattering phase shifts and amplitudes of heavy 
and light hole wavefunctions. The Green's matrix of the Schrodinger equation, the Lippmann-Schwinger equation and 
the Born approximation are discussed in detail in section IV. The object of Section V is the study of the non-trivial 
behavior of the wavefunctions near the origin, the knowledge of which is crucial for successful numerical computation 
of the solutions of the abovementioned differential equations. In section VI we give examples based on Yukawa-type 
potentials to illustrate our calculations before concluding. 



II. SPHERICAL MODEL HAMILTONIAN IN THE STRONG COUPLING REGIME 

Luttinger showed that in the strong coupling regime the holes can be seen as quasi-particles of spin J = 3/2 [26| . 
Baldereschi and Lipari treated the problem of acceptor states in the effective mass approximation in order to separate 
the symmetrically spherical and cubic terms in the Hamiltonian They showed that the spherical Hamiltonian, 
7i s , may be written as: 

W.(r,p) = P 2 - „ ■ + V(r), (1) 

Zttiq lemo 

where mo is the free electron mass, // = (472 + 673)7571, and 71, 72 and 73 are the Luttinger parameters j2(|. The first 
term of TL S in Eq. Q is the kinetic energy, the second term represents the spherical part of the spin-orbit interaction, 
and the last one is the perturbation potential created by the defect that is assumed to depend only on the distance r. 
The diagonalization of the unperturbed Hamiltonian, 7i® (p), where p = ?ik, yields two parabolic bands characterized 
by the effective mass coefficients 1/(1 + /i)7i = r 2 /7i for the light holes and 1/(1 — = ^m/71 for the heavy holes. 
Denoting E(k) — fi, 2 7ik 2 /2mo the hole energy which is a constant of motion, the wavevectors of the light and heavy 
holes are respectively k + = r p k and k~ = r m k. 

The Hamiltonian Ti. s commutes with the total angular momentum F = L + J, where L is the angular momentum. 
Therefore the scattering eigenstates may be written as 1 3/2, F, M), where F(F + 1) and M are the eigenvalues of F 2 
and F z respectively. 
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TABLE I: Values of cos a and sin a as function of £ and L. 



f — — 1L — F+i L > 2 


f — 1 L — F — i L>1 


cosa- 


COSQ "2T+T 


„ ina _ VHL + l)(L-l) 
sma 2L + 1 


„ ina _ >/3L(L + 2) 
sma 2L + 1 



For a given value of F (F > 3/2) there are 4 corresponding values of L: L = F — 3/2, L = F — 1/2, L = F + 1/2 
and L = -F + 3/2. Moreover since the parity II = (— l) L is conserved, only the states \L,3/2, F, M) satisfying 
AL = 0, ±2 are coupled. In the sub-space \F, M), the Hamiltonian matrix Ti s contains two 2x2 blocks of well 
defined parity. The first block is associated to the states {\L = F — 3/2), \L = F + 1/2}}, and the second to the states 
{\L = F — 1/2), \L = F + 3/2}}. From now on we denote \L— 1) and \L + 1) the two coupled states of parity II and we 
shall use the parameter £ = ±1 introduced by Ralph [ll| to distinguish the value of L corresponding to each block: 
if £ = 1, L = F - 1/2 and if £ = -1, L = F+ 1/2. 

The two radial wavefunctions /i±i(r) associated with the scattering states \F, £,M) satisfy a system of coupled 
radial Schrodinger equations, which reads |l2j : 



Hl-ix-i F[l-i,l+i 
Hl+i,l-i Hl +1 ,l+i 



Jl+i 



E 



Jl-i 
Il+i 



where the four matrix elements -Hl±i,l±i ar e second order differential operators of the variable r. 
potential V(r) appears only on the diagonal of the Hamiltonian matrix. 



(2) 

The perturbation 



III. DERIVATION OF THE PHASE EQUATION 

Ralph has demonstrated that in a region where the perturbating potential has no effect, the two coupled radial 
wavefunctions fL±i{f) may be characterized by four constant parameters: two amplitudes a ± and two phase shifts 
5^ [TT| . In this section we show how we obtain those four parameters as the limits, when r — > oo, of the functions 
a ± (r) and associated to a spherically symmetric potential truncated at distance r. 

A. From 2x2 second order differential system to 4x4 first order differential system 

To eliminate the factor in the solutions of the unperturbed system and in their asymptotic expressions, we 
define the functions u_L±i(r) = rfi,±i(r). Introducing the dimcnsionlcss variable x = kr and the reduced potential 
v(x) — V/E, and factorizing the differential operators Hl±i,l±i, we find that ul±i(%) satisfy the following system: 



(1 + /icosa) 
— fi sin a 



_d_ 

dx 



d , L 
L + l 



— /ism a 
(1 — fx cos a) 



d , L 



d _ L 
Hi x 



\d L+l 


A 


[ d , L + l] 


[dx x 




[dx + x 




= (v(x) - 1) 



UL-l(x) 

u L+ i(x) 



(3) 

where the angle a introduced by Gel'mont and D'yakonov |2j| is defined for each value of F and £ in Table 1. 

To lower the order of this differential system and to obtain the free solutions (v = 0) in a simpler way, we define 
the functions G (x) as follows: 



\=[B\-'[tt\- 



-1 
1 



d _ L 
die x 



UL-l(x) 



d 
die 



L + l 



u L+ i(x) 



(4) 



where [f2] is the 2x2 rotation matrix of angle a/2 and [R] the diagonal matrix whose elements are r p and 
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Introducing the vectors u = y " L 1 J and G = ( ^_ J , the system in Eq. take the form of a hrst order 
differential system perturbed by a 4 x 4 matrix [V(x)]: 

( " N [A(x)](i)-[V(x)}(%), (5) 



dx \ G J \ G J \ G 

which is a very convenient approach to study the solutions near the origin. In Eq. J5J, L4(x)] an d [V(x)] are defined 
as follows: 



[A(x)}= 



i o \ [nr 1 ( -l o 



(6) 



and 



P\ \0] 

where [O] is a 2x2 matrix whose elements are all zero. 

B. The phase equation 

When the potential v(x) is zero, we find that the vector G satisfies a differential equation of the Schrodinger type: 



^ X L(L + lj - r „, 2 

d.i: 



2 G — 2 — G + [R] G = 0. (8) 



The above equation has four independent solutions of which two are regular solutions G^ r and two irregular G^;. 
The vectors and , respectively associated to Gq and Gq , may be readily obtained from Eq. JSJ with v = 0. It 

is convenient to put together the four free solutions, regular (r) and irregular (i) ( ^ ) into a 4x4 matrix [W(x)] 

V G /0,r/i 

whose elements are given in Appendix A. The regular and irregular parts of this matrix may be expressed in quite a 
simple fashion with modified spherical Bessel ji and Neumann hi functions respectively |Tfij . The general solution for 
v(x) = can thus be written: 



= [W(x)} Go (9) 



where Go is a four-component constant vector. 

When the potential v(x) is non-zero, we apply the Lagrange method of the variation of constants to find the 
solutions: 



= [W(x)}C(x). (10) 



Since [M^(x)] is solution of the unperturbed system, i.e. d [W] /dx — [A] [W], we find that the unknown vector C(x) 
satisfies the following differential equation: 



dG 
dx 



[w]~ 1 [v][w]c = -[wr 1 [v]^y (11) 
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In the framework of the variable phase method [T^ the vector C(x) is searched in the form: 




C= h (12) 



where, for ease of notation, we have omitted the quantum numbers (-F, £, L). 

Inserting Eq. (fX"Sf) into Eq. (|TT]l and performing the algebra yield the following system of differential equations: 

A ( ^) A = _ v(x) ( cos ^ " L - l(rpX) sin ^ ^ ( ^ (m 

dx V S+ W J { cosf ^_!(r p x) sinf j i+1 (r p2; ) J V "^W ^ ' 

^ /-sinf n L _ l( r m .) cos § ft L+1 (r m x) \ /„ w( ^ 

d * V 7 ^ - sin § 5x-i(r m x) cos § j L+1 (r m x) J V 7 " 

The variable phase shifts d ± (x) and amplitudes a ± (x) are now introduced in the same fashion as in problems 
without coupling |l5j : 

(x) = (x) cos 6 ± (x) , (14a) 



s ± (x) = a ± (x)smS ± {x). (14b) 

Inserting the above expressions of and in the differential systems Eqs. (|13a|l and l)13bjl . we obtain the 
generalized non- linear coupled equations for the phases and amplitudes of the holes: 



-r v D + +r m a D c /a + 
dx \ a+ | - V[X> | -r p a+SD+ +r m a-SD+> 



(15) 



—r m a SD + r p a + SD c ,+ 



where the functions , 5D ± , D c and SDp~ £ (e = ±) are given in Appendix B. 

Since the wavefunctions u_l±i must be regular at the origin, the above differential system must be solved with 
<5 ± (0) = as initial conditions. This choice is justified by the fact that if the potential is truncated at the origin, the 
perturbation is removed from the problem and hence the shifts of the phases of the wavefunctions are zero. However, 
the behavior of the functions a + and a~ near the origin is not imposed. There are thus two degrees of freedom that 
allows one to generate two independent solutions. This point is discussed further in Section V. 

C. Asymptotic expressions of the wavefunctions 

We assume that the perturbation potential V(r) has a finite range rg (the Coulomb potential is therefore excluded). 
For x > kro and x S> L/r v (y = p,m) the functions a ± (x) and S ± (x) are practically constant: 

c^(x) w a ± (oo) cosi5 ± (oo) = a cosfJ* (16a) 



s ± (x) « a ± (oo) sin(5 ± (cxD) = a ± sin<5 ± , (16b) 
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and the functions ji(r v x) and hi(r u x) can be taken as circular functions: lim^—,,^ ji(r v x) — sin [r u x — Itt/2) and 
linia^oo hi{r v x) = — cos (r v x — Itt/2). Using Eq. (|1(J|) and formulas of Appendix A, we find the asymptotic expressions 
of the wave functions wl±i : 

OL OL 

u l-i(x) ~ r p a+ cos cos (^p^ — Ltt/2 + J + ) — r m a~ sin — cos (r m x — Lir/2 + S~) (17a) 

OL OL 

« — r p a + sin — cos (r p a; — Ltt/2 + <5 + ) — r m a~ cos — cos (r m ai — Ltt/2 + <5~) . (17b) 

Comparing these expressions to those derived in Ref. we see that the variable phase method is suitable to 
compute Ralph's parameters. 

The functions ui and Uh defined as linear combinations of ul-i and ul+i'. 

u f\)=[n]( UL -f\) (is) 
uh(x) j y u L+ i(x) j v 1 

exhibit a sinusoidal behavior in the asymptotic region, i.e. x — > oo: 

u/{x) w r p a + cos (r p x — L7r/2 + 5 + ) (19a) 
Uh(x) ~ — r m a~ cos (r m x — Ltt/2 + <5~) (19b) 

Given the definitions of r m and r p in Sec. II, we see that ug and Uh are the light and heavy hole radial wavefunctions 
respectively. 

IV. BORN APPROXIMATION 

A. Green's matrix of the Schrodinger equation 

Inserting the expressions of c^(x) and s ± (x) obtained from the integration of Eqs. I|13a|) and l|13b|) . in Eq. |JTDJ, 
yields the Lippmann-Schwingcr integral equation in terms of partial waves [29j, satisfied by the two components ul±i 
of the regular solution: 

u L -x{x) \ _( <_x(a-) \ [ x ' ( f G°l-i.l-i(x,x') Gl-LL+^X.X 1 ) \ f ul-^x 1 ) \ , , 
u L+1 (x) ) ~ { ul +1 (x) ) + J V ^ X >\ gl +l . L ^{x,x') Ql +W {x,x>) ) \ u L+1 {x>) ) ' ^> 

where the four matrix elements Gl±i l±i( x > x ') are gi ven m Appendix C. The functions vP L:ill are the components of 
a regular solution of the Schrodinger equation of a free particle. They are obtained by taking = and = in 
Eq. HlOfl . where and Oq are two constants that appear in the integration of the differential system, Eqs. (|13a|) and 
(|13bjl (see Appendix A). 

It is also possible to show that the matrix \Q°(x, x')] whose elements are Gl±i l±i( x > x ') ^ s a Green's matrix of the 
hole unperturbed Hamiltonian, [-ff ] ■ Indeed, if one writes the Schrodinger equation, Eq. 0, as: u = r m r p v ( x )'u, 
one checks that the matrix \Q°~\ satisfies the following equation: 

[H°] [G°] = rlrl 5{x-x')[\\y (21) 

where the Dirac (^-function comes from the differentiation of the step function @{x — x') that appears in the expression 
of the four matrix elements of [G°] (see Appendix B). 
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B. Born approximation 

The partial waves may be obtained in the Born approximation if the functions u < [ i ± 1 (x') are substituted to the 
functions Ul±x(x') in the right hand side of the Lippmann-Schwingcr equation, Eq. I|2()|) . 

It is possible to find approximate expressions for the functions S ziz (x) and a^(x), by comparing uf^^(x) to the general 
expression of Ul±i(x) given by Eq. (|10|l . More precisely the identification of the terms in fiL±i(r p x), fiL±i(r m x), 
jL±i( r px) and jL±i{ r mx), gives the functions S + (x), 5~(x), a + (x) and a~(x) respectively. Noting that 5 ± <§; 1 and 
keeping only terms of the first order in the potential strength yield: 



S+{x) 



v(x')(cos 2 I jl.iOyr') +sin 2 | f L+1 (r p x')^dx' 



v{x')F{x')dx' 



(22a) 



r(i)«-r m f\(x')(sm 2 ^ f^i^x') + cos 2 ^ f L+1 (r m x'))dx' - r p ^ f v(x')F(x')dx' (22b) 
Jo v 2 1 ' a-n Jo 



where 



F(x) =sin^ cos^ x [j L+1 (r p x)j L+1 (r m x) 



h-l{rpX)j L -i{r m x)] 



(23) 



and 



a + (x) 



1-7Y 







a a 
a r m sin — cos — 



(x')(cos 2 ^ h L _i(r m x')] L ~i{r p x') + sin 2 ^ n L+ i(r p x')j L+ i(r p a;')^dx' 

w(a; / )("L+i( r P a; ')ii+i( r ms') - ni-i(r p a;')j L _i(r m x'))da;' (24a) 



1 - r m y w(a;')^sin 2 - n L ^i(r m x')j L -i(r m x') + cos 2 - h L+ i(r m x')j L+ i(r I[L x r )jdx' 
(x (y. f x / 

- a(j> p sin- cos- y w(xO^ L+ i(r m a;')iL+i(? , p^ / ) ^ ^i-iC^m^OiL-^^p^OJda;', (24b) 

Since n; (x) » -(2/ - 1)!! a;"' and j;(ir) « x' +1 /(2Z + 1)!! near the origin 0, the integrals in Eqs. l(2^al) . (2201, 
(|24a|l and (|24b|l . converge for potentials that vary like x~ a with a < 2 near x = 0. The presence of terms in (aj ) _1 
and (dp - ) -1 in i5 + and <5~ is not a problem since it is always possible to generate two independent solutions taking 
ciq =/= and ^= 0. We note that since the phase shifts depend on the ratio /a^ , which is arbitrary, there is no 
correlation between their sign and the attractive or repulsive nature of the potential. Therefore the meaning of the 
phase shifts 5^ is not as clear as it is for the phase shift Si of the partial waves when there is no coupling. 



a (x) w a 



V. BEHAVIOR OF THE WAVEFUNCTIONS NEAR THE ORIGIN 

The right hand side of Eq. (|15() satisfied by the scattering amplitudes a ± (a:) and phase shifts S ± (x), contains 
Neumann's functions hi(x) that diverge when x — > 0. In some cases the reduced potential v{x) may also diverge. 
Therefore, the computation of the solutions of Eq. (|15|l requires the knowledge of the behavior of either the functions 
a ± (x) and 5 ± (x) or the functions c ± (x) and s ± (x). Equations l|13afl and l|13bjl show that this problem is equivalent 
to finding the approximate expressions of the wavefunctions ul±i(x) valid near the origin. 

One method is to solve the integral Lippmann-Schwinger equation, Eq. (|2U|) . by successive iterations. The most 
difficult case that may be treated in the present work is that of a potential v(x) diverging in x — as a Coulomb 
potential, i.e. v(x) w Vqx^ 1 when x — > 0. Singular potentials behaving like x~ m with m > 2, near the origin are not 
considered in the present work. After the first iteration, we see that all the integrals appearing in u^iiC^) converge 
in x = 0. However, in the general case, it is impossible to numerically solve the differential system Eq. (|15J) by 
substituting u^^x) to ul±i(x) near the origin. It is necessary to go to order 2 in the potential strength. After the 
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second iteration, we note that one integral of u^^x) diverges logarithmically like Vq x' 1 dx' unless and a 
satisfy: 

4 cos I r£ +1 = ao sin | r^ +l (25) 

If the above equality is not satisfied, it is not possible to solve the differential system Eq. I) 1511 by imposing the values 
of the amplitudes at the origin, a ± (0). 

The general behavior of the regular solutions of the Schrodinger equation near the origin was studied by Newton 
in the case of a tensorial scattering potential, represented by a 2 x 2 matrix [3(J • For this kind of coupling the Green's 
matrix is diagonal. Using the Lippmann-Schwinger equation, it was found that one solution can be expanded as 
power series whereas the other exhibits an amplitude function that diverges logarithmically from the first iteration. 
Newton showed that one could circumvent this difficulty by modifying the inhomogeneity of the Lippmann-Schwinger 



equation, in our case o _1 I , in a convenient way so that the problematic terms cancel each other out |30l l31j. 
V u l+i J 

According to Newton this method may always be applied in a general case, but the way to modify the integral equation 
may be more complicated either if the difference between the values of the two coupled angular momenta is greater 
than 2 or if the coupling is of a different nature. 

The case of the scattering of holes is indeed very different: the potential is scalar, the Green's matrix is not diagonal 
and diverging terms appear only at the second iteration. This renders any use of Newton's method very difficult. 
For this reason we have adopted another method, which consists of taking the first order differential system, Eq. JSJ, 

satisfied by the vector ( ^ defined in section III. A, as a starting point. The problem of logarithmic divergence that 

we found in the present problem, has actually its origin in the presence of a a; -1 term, called singularity of the first 
kind, in the expansion of the matrix [j4(a;)] — [V^rr)] near the origin. One result of the theory of differential systems 
is that the general regular solution generated by two arbitrary constants (a, b) may be expanded as follows |32j |: 



(26) 



where i and j are integers that a priori satisfy: < i < 3 and j > 0; r\ is an eigenvalue of the constant matrix 
L4o] that is the first term of the series expansion of the matrix [^4(x)] defined in Eq. JfjJ): [^4(x)] = J2 n >o i-^n] % n ■ 

The four-component constant vectors Ci.j(a, b) are determined recursively from the series expansion of [^4(a;)]. In the 
present work, the problem is reduced to the diagonalization and inversion of 4x 4 matrices. 

The relationships between the constants (a, b) and (o^jOq) can be easily obtained from the comparison of the 
lowest order terms of the series expansions of M( a ,b) and u° L±1 {x). We found that that there exists one solution that 
may be expanded as a power series, generated by the couple (a — 0,6) |33|. At the lowest order the wavefunctions 
behave like: ul±i{x) » x l+2 . One can check that for this solution the equality in Eq. (|25|l is satisfied. If a ^ 0, the 
expansion contains a power series plus an additional term originating from the product of Vq log(x) and the previous 
solution obtained with (a = 0, b = 1). The log(x) term, known in the theory of Fuschian differential equations 
appears only at the second order for potential behaving like a Coulomb potential near the origin, i.e. v(x) w Vqx~ 1 
when x — > 0. The lowest term of the power series varies like x L for ul-i(x) and Vqx l+1 for ul+i{x). Note that the 
functions wl±i are regular at x — for any couple (a, b), despite the presence of a log(:c) term in the expansion. 

The functions c ± (x) and s ± {x), and hence a ± (x) and ^ ± (x), are deduced from ul±i(x) and G^t^x) solving Eq. 11UI) 
(we found that det[IU] = r m r p ). For any solution generated with o^0, the constants Oq do not correspond to the 
initial values of the amplitudes a ± (x) since these functions diverge like log(x). They only coincide with the asymptotic 
values a ± (oo) if they can be calculated in the Born approximation for a very weak potential. For the solution that 
can be expanded as a power series, the constants and a,Q can be identified to the initial values of the amplitudes 
a + (0) and a~(0). In the Born approximation, the equalities a ± (oo) ~ are, of course, still valid. 

The two independent solutions that we have just described exhibit the same behavior that we observed on the first 
iterations of the Lippmann-Schwinger equation at the beginning of this section. The power series appearing in the 

expansion of ( ^ have a finite convergence radius £conv that limits its use to a restricted interval; moreover x conv 

is in general smaller than the point where the functions (x) and 5 ± (x) cease to vary and which is directly linked to 
the range of the potential v(x). To obtain a ± (oo) and 5 ± (oo), we thus solved the differential system Eq. I|15[l using 
the Runge-Kutta method from a value of x, x = x m < x conv . 
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VI. NUMERICAL CALCULATIONS 



Numerical calculations were performed to solve Eqs. (|15fl for the scattering phase shifts, amplitudes and radial 
wavefunctions. As examples we considered the scattering of the wave in the state \F = 3/2, L — 1), in Yukawa-type 
potentials frequently used to model charged point defects induced perturbation. 



A. The Yukawa potential 



The Yukawa potential we chose for the present example reads: 



v(x) = , 27 

ka-Q x 

where Z is the charge number, ob the Bohr radius, and A s the screening length. 

The potential above corresponds to the linearized Thomas- Fermi (LTF) description of the charge density devoted to 
the screening of a sing le ion. The LTF theory neglects the many-body effects contained in the Lindhard or Hubbard- 
Sham formalisms 36] . The assumption of independent scattering by the various ions is valid only as long as the 
average defect separation is large compared to both the hole de Broglie wavelength and the screening length so that 
neighboring ion potentials do not overlap significantly |37| . 

In the Born approximation, the screening length A s is equal to the Dingle-Mansfield value: A~ Born = 
q 2 /eot r k^T dp/AE-p [23], where p is the hole density, Ef the Fermi energy and e r the relative static dielectric constant. 
When the Born approximation is invalid, the energy dependence of the screening charge must be taken into account. 
In this case the screening length A s is evaluated numerically and self-consistently using the generalized Friedel sum 
rule ,3^ which insures that a given ion is fully screened at large distance. 

Numerical calculations were performed taking the following values: Z — — 1, A 8 = 20ob with = A^t^t^^xj ra^q 2 
and /z = 0.481 for the effective mass parameter of Si ^jj. 



1. Solutions without logarithmic divergence term 



Scattering amplitudes and phase shifts of the solution without logarithmic divergence term are depicted in Figs.^ 
and|2 for various values of the energy of the incoming hole, E(k) = (fcfle) 2 -Eo, where Eq = ft 2 7i/2moa B is the Rydberg 
energy. The curves were obtained by solving the differential system of Eq. 1)15(1 using the Runge-Kutta method, with 
a~(Q) = ao = ( 2L + 3 ) !! / r m +1 cos(a/2) and a+(0) = = (2L + 3)!!/r£ +1 sin(a/2) taken as initial conditions (see 
Eq. 12510 . This choice corresponds to the following values: <Jq = 5.505 and = 15.708 for the \F,L) state that we 
consider. 




x = kr x = kr 



FIG. 1: Scattering amplitudes a (x) and a + (x) as function of the scaled distance x = kr, for various values of kaB- 



Note that for fcae > 100 the phase shifts become very small and the amplitudes vary very little from their initial 
values. This is the regime where the Born approximation is valid and may be applied. Note that some of the phase 
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x = kr x = kr 

FIG. 2: Scattering phase shifts 5~ (x) and S + (x) as functions of the scaled distance x = kr, for various values of has- 

shifts on Fig. [2] have not quite reached their asmptotic value for x = 30. This is of no importance here since these 
quantities are not further exploited. 

The radial wavefunctions urj(x) and U2(x) as well as ug(x) and Uh(x) of the solution without the logarithmic 
divergence term in Fig. [3] were obtained from the scattering amplitudes and phase shifts discussed above. 




10 20 30 40 10 20 30 40 

x = kr x = kr 



FIG. 3: Radial wavefunctions uo(x) and U2(x) (left) and light and heavy holes wavefunctions ue(x) and Uh (right) computed 
for fcdB = 0.15. 

One can check on Fig. [31 (right panel) that the light and heavy holes wavefunctions ue(x) and u^x) already start 
varying in a quasi-sinusoidal fashion when x approaches x = 40. 

The influence of the orbital momentum L and the charge number Z of the defect on the phase shifts <5 ± (solution 
without logarithmic diverging term) is depicted on Fig. 0] Only the states £ = 1 for which L = F — 1/2 were 
considered. We took the following values for the parameters: Z = — 1, A s = 20ae and ka-e. = 0.5 for the study of the 
influence of L on and L = 1, A s = 15ae and fcae = 1 for the study of the influence of Z on 5 ± . 

One striking feature is that the phase shifts S + (x) and 8~(x) behave very differently. The function 6 + (x) varies 
monotonically, is positive for attractive potentials and negative for repulsive potentials. Moreover its value decreases 
with increasing orbital momentum L. The function 5 + (x) has all the properties of the phase shifts of a scattering 
problem without coupling |15|. The phase shift S~(x) exhibits essentially two atypical behaviors: it is negative for 
attractive potentials (Z — —1 and Z — —2) and its behavior as a function of L is "normal" only for small values of x. 
Indeed, outside the small x region, 8~{x) strongly oscillates before adopting a new trend: the absolute value of 5~{x) 
increases with increasing orbital momentum L. 

This unexpected behavior of 5~(x) can be explained in the light of the study of the coupling terms a~D c /a + 
appearing in the expression of the derivative of 5 + (x), d<5 + (a:)/da;, and a~D c /a + appearing in the expression of 
d<5~(cc)/d:r (see Eq. O). The study of the ratio a + (x)/a~(x) (not given here) shows clearly that, for a given value 
of L that we take here to be L — 4, it decreases for x < 5 before increasing suddenly, oscillating in the vicinity of 
x w 10 and finally converge towards 20. This coupling term has therefore a strong influence on the behavior of S~(x), 
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" 6 5 10 15 20 25 " u 10 20 30 40 50 60 
x = kr x = kr 

FIG. 4: Scattering phase shifts S~(x) (full line) and 8 + (x) (dashed line) as functions of x for various values of orbital momentum 
L (left panel) and charge number Z (right panel). 

while it does not affect 5 + (x) significantly. 

2. Solutions including logarithmic divergence term 

Another solution whose amplitudes presents a logarithmic divergence, was generated. We chose the following values 
for the integration constants: = — 1 and = 1, and kept the previous parameters: A s = 20 a^, Z = — 1 and 
L = 1. The logarithmic divergence of the amplitudes near x = appears very well on Fig. For large values of the 
energy of the hole, E(k) (fcae > 100 for a~ and fcae > 10 for a + ), the asymptotic values, a ± (oo) and the parameters 
are practically the same. 

o 



-0.2 
-0.4 
_-0.6 

X 

-0.8 



-1 



-1.2 



5 10 15 20 5 10 15 20 
x = kr x = kr 

FIG. 5: Scattering amplitudes as functions of x considering the logarithmic divergence term, for various values of fcae- 

The radial wavefunctions uq(x), U2(x), ug(x) and ith(x), obtained for fcae = 0.15 are depicted in Fig. EI 
One can check that the logarithmic divergence does not compromise the regularity of the radial wavefunctions near 
x = and that ui{x) and Uh(x) start exhibiting sinusoidal oscillations from x « 30. 

B. Modified Yukawa potential 

The last case treated here is that of the scattering of the wave \F = 3/2, L — 1) by a repulsive Yukawa potential 
(Z = 1, A s = 200 ae) modified by an attractive square core (V — — 1.6E for r < n = 10 ae). This modified Yukawa 
potential has a positive maximum V^nax ~ 0.19 Eq. 
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x = kr x = kr 



FIG. 6: Wavefunctions uo, 112, ue and tth as functions of x. 



Because of its rapid variations, this type of potential is not rigorously compatible with the use of the effective mass 
approximation. It is nonetheless frequently used to model the chemical shift that is observed on the bound states 
of doping impurities |4l| . It is not in the scope of the present work to study the effect of the core potential on the 
phase shifts 6^, but it is of interest to put in evidence the influence of the energy E(k) on the radial wavefunctions 
Ul±i(x) in the case of a non-monotonical potential. Since the perturbation potential is regular in x — (Vb = 0), 
the amplitudes do not exhibit a logarithmic divergence and the two independent solutions, (a) and (b), can be 
expanded as a power series. 

The behavior of the scattered waves is different if the energy of the incoming hole, E(k), is greater or smaller than 
Vmax- We chose two values for the hole energy: E(k) — 0.64 V max and E(k) = 4 V max , that corresponds to fcae = 0.349 
and fcae = 0.872 respectively. The reduced potentials and hole energies are depicted in Fig. J7J). 



1 1 1 1 r 



(b) 



0- 



-12 - 



_J I I I I I L 

5 10 15 20 5 10 15 20 

x = kr x= kr 



FIG. 7: Modified Yukawa potentials as functions of x. The dotted lines indicate the location of energy of the hole with respect 
to the potential barriers in full lines. Two cases are considered: v(x) = V/E(k) for E(k) — 0.64 Vmax on figure (a) and 
E(k) = 4 V max on figure (b). 



The radial wavefunctions uq(x) and 1/2(2;) are depicted on Figs. 8 and 9. They were obtained by taking the couple 
(a-(0) = a = cos^r /r£ +1 = 0.367, a+(0) = a+ = sin ^ /r£ +1 = 1.047) for the solution (a) and (or(0) = = 
— 1, a + (0) = = 1) for the solution (b) as initial conditions. 

For E(k) < Vmax the vanishing wave that goes through the classically forbidden region, [3.5; 5.5], by tunnel effect 
can be seen on the radial wavefunctions uq(x) and U2(x) of the solution (b), in Fig. EI The high value of r\ allows 
one to clearly see the rapid oscillations in the region where the potential is attractive. 

For E(k) > Vnax (Fig- the wavefunctions exhibit an oscillatory behavior for all the values of x and their 
amplitudes are of the same order of magnitude both in the well region and the barrier. 
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5 10 15 20 
x = kr 



FIG. 8: Radial wavefunctions uq(x) and U2(x) of the solutions (a) and (b) for E(k) = 0.64 V m 

0-4 3 r 

0.3 - " 




10 15 20 
x = kr 



FIG. 9: Radial wavefunctions uo(x) and U2(x) of the solutions (a) and (b) for E(k) = 4 Vm 



VII. SUMMARY, DISCUSSION AND CONCLUSION 



In this article we addressed the problem of realistic calculations of scattering states of holes in mixed valence bands 
in the strong spin-orbit coupling regime. From the early stage of the present work it appeared that the generalization 
of Ralph's results to finite range potentials, such as those created by ionized defects, was indeed specially involved. In 
particular, the behavior of the two radial wavefunctions ul±i{x) near the origin required special attention. We found 
that the variable phase method provided us with a very convenient framework that facilitated such calculations. 

Starting from the two coupled Schrodinger equations satisfied by the radial wavefunctions, and introducing auxiliary 
functions, we first derived the system of four coupled non- linear differential equations whose solutions are the two 
phases and two amplitudes characterizing the scattering states with well defined total angular momentum and parity 
\F,L). This systems forms the generalized phase equations [l5l . Ralph's parameters were obtained by integrating 
these equations over an interval larger than the range of the perturbation potential. Then we obtained an expression 
of the Green's matrix of the unperturbed Schrodinger equation as well as the Lippmann-Schwinger integral equation 
which yields an iterative scheme to evaluate the radial wavefunctions. Considering only the first iteration, we found 
the "Born approximation" for the amplitudes and phases of the scattering states. 

Owing to the non-diagonal form of the Green's matrix, the expansions of the solutions near the origin (x ~ 0) that 
are necessary for successful numerical integration of the generalized pha se equation, have been constructed by making 
use of a technique only applicable to first order differential systems [33 • We showed that for Coulomb- like potentials, 
there exists logarithmically diverging terms in the amplitudes near the origin, appearing only at the second order 
in the potential strength. This specific result stems from the coupling of the holes wavefunctions: for the tensorial 
couplings that we found in the litterature 

HHUm, the log(x) terms are indeed generated from the first iteration 
of the Born series. This fact might explain why we could not find published results on the problem we treat in the 
present paper since Ralph's article [Tl| . 
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The generalized phase equations were solved for screened Coulomb potentials of the Yukawa type and various values 
of angular momentum, incoming hole energy and defect charge state. The oscillatory but non-sinusoidal behavior of 
the radial wavefunctions ul±x(x) in the asymptotic region is coherent with the predictions of Ralph's theory [Tlj . 
By linear combination of those functions, we constructed the heavy and light hole wavefunctions which exhibit a 
sinusoidal asymptotic behavior. This result shows that the phase shifts, solutions of the generalized phase equations, 
are those of the heavy and light holes. However their meaning is not as clear as for the phase shifts Si of the partial 
waves in the absence of coupling since their signs do not only depend on the charge of the ionized defect but also on 
the choice of two arbitrary integration constants. 

Calculations of Hall and drift mobilities require proper computation of the energy average of the various relaxation 
time functions 0, 0, El • A study of the influence of the temperature and concentration of doping atoms on those 
physical quantities is practically achievable only on the condition that the variations of the relaxation times , over 
a range of incoming hole energy E and screening length A s that is accessible to experiments, could be reproduced by 
approximate analytic forms analogous to those of Meyer and Bartoli 0. A vast amount of purely numerical work 
remains to be done to construct those approximate analytic expressions because it may appear necessary to consider 
a large number of scattering states \F,L) for some values of (E, A s ). This work is beyond the scope of the present 
article. 

The strong spin-orbit coupling regime is well adapted to semiconductors like Ge and GaAs, where the energy of 
spin-orbit interaction is about Ago ks 300 meV. The study of the scattering states in the weak coupling regime, i.e. 
Aso ~~ * 0, can be seen as the other limit case that can be tackled before addressing the more involved scattering 
problems in Si for which Aso = 40 mcV, i.e. well in the intermediate coupling regime. The variable phase method 
can be used again since the two radial wavefunctions u^±i are solutions of a differential system of the second order 
analogous to the one in the strong coupling regime [H| . 



APPENDIX A: FREE SOLUTION 



The 4x4 matrix [W(a;)] whose columns are the four linearly independent free solutions (v = 0) of the differential 
system Eq. (JSJ, reads: 



[W} = 



'-o.i 



'O.i 



f+ c+ r<~ c 

b 0,r °0,i °0,- Lt ' 



O.i 



(Al) 



where 



r+ _ ( h(r P x) 

q+ = ( n L (r p x) 



lUldG °>*-\ fr(r m x) 



and G ^-\n L (r m x) /' 



(A2a) 
(A2b) 



and 



rpCOSTj z L ^x(r p x) 
r p shiTj z L+ i(r p x) 



(A3a) 



-r m sin^ z L _i(r m x) 
cos y z L+1 (r m x) 



(A3b) 



where ji and hi are the modified spherical Bessel and Neumann functions respectively |l5j and zi = ji for the regular 
solutions and ii — hi for the irregular solutions. 
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Taking C = 







in Eq. yields the free regular radial wavefunctions that appear in the Lippmann-Schwinger 



V o / 

equation, Eq. lj2U|) : 



u l-i( x ) = r P cos \ a>oJL-i(+) - r m sin ^ a j L -i(~) 



(A4a) 



u° L+1 (x) =r p sin- a h j L+1 {+) + r m cos - a j L+1 (-) 



(A4b) 



APPENDIX B: FUNCTIONS D ± , SD ± , D c AND SD E C 



Below are given the various functions that appear in the generalized equations for the phases and amplitudes, 
Eq. fT5|l: 



^=cos 2 f [^ 1 ] 2 +sin 2 ^ [Dt ±1 ]\ 



(Bl) 



SD ± = cos 2 | S± T1 £± T1 + sin 2 St ±1 D$ ±1 , 



a 



(B2) 



D C = D+_ 1 D-_ 1 -D+ +1 DZ +1 , 



(B3) 



with e = ±, 



and 



SDl>-e = sin | cos | {Si^DZ^ - Si+^i) 



£>± = ji(±) - sin(«5 ± )n i (±), 



Sf = cos( ( y ± )n z (±) + sin(5 ± )j i (±). 



(B4) 



(B5) 



(B6) 



For ease of notation in Eqs. (|B5|) and (|B6(I . the arguments r m x and r p x of the modified spherical Bessel and 
Neumann functions are denoted — and + respectively. 



APPENDIX C: GREEN'S MATRIX ELEMENTS 

The expressions of the four matrix elements of the Green's matrix [Q°~\ defined in Eq. (|20|l : 



Gl-1,L-i( x > x ') = r P cos2 -^ 9L-i(r p x,r p x') +r m sin 2 - gr L _i(r m a;, r m x') 



&l+i,L+i( x > x ') = r p sin2 f 9L+i(r p x, r p x') + r m cos 2 ^ g L+ i{r m x, r m x') 

Gl-i, L +i( x , x ') = sin ^ cos^ r p (jj L+ i(r p x')n L ^i(r p x) - j L -i(r p x)n L+ i(r p x'f) 

- r m (jL+i(r m x ')nL-i(r m x ) - h-i{ r mx)n L+1 {r m x')^ Q(x - x') 



(CI) 
(C2) 

(C3) 



- r m (j L -i(r m x')n L+1 (r m x) - jL+i(r m x)h L - 1 (r m x')j Q(x ~ x'), 
where @(x) is the step function and 



gi(x,x')= ji{x')hi(x) - ji(x)ni(x') xQ(x-x'). 
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